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ABSTRAGT 

Interpreting the spectra of brown dwarfs is key to determining the fundamental physical and chem¬ 
ical processes occurring in their atmospheres. Powerful Bayesian atmospheric retrieval tools have 
recently been applied to both exoplanet and brown dwarf spectra to tease out the thermal structures 
and molecular abundances to understand those processes. In this manuscript we develop a signifi¬ 
cantly upgraded retrieval method and apply it to the SpeX spectral library data of two benchmark late 
T-dwarfs, G1570D and HD3651B, to establish the validity of our upgraded forward model parameter¬ 
ization and Bayesian estimator. Our retrieved metallicities, gravities, and effective temperatures are 
consistent with the metallicity and presumed ages of the systems. We add the carbon-to-oxygen ratio 
as a new dimension to benchmark systems and find good agreement between carbon-to-oxygen ratios 
derived in the brown dwarfs and the host stars. Furthermore, we have for the first time unambiguously 
determined the presence of ammonia in the low-resolution spectra of these two late T-dwarfs. We 
also show that the retrieved results are not significantly impacted by the possible presence of clouds, 
though some quantities are significantly impacted by uncertainties in photometry. This investigation 
represents a watershed study in establishing the utility of atmospheric retrieval approaches on brown 
dwarf spectra. 


1. INTRODUCTION 

The spectrum of a brown dwarf opens a series of win¬ 
dows into the depths of its atmosphere, revealing its com¬ 
position and thermal structure. Differing wavelengths 
peer to a diversity of depths and are influenced by varying 
atmospheric species. Teasing out the chemical composi¬ 
tion and atmospheric structure of brown dwarfs by syn¬ 
thesizing the information emerging from each window is 
key to understanding the processes acting in their atmo¬ 
spheres and possibly their evolutionary histories. Histor¬ 
ically, empirical analysis techniques, such as color mag¬ 
nitude diagrams, the appearance of spectroscopic fea¬ 
tures, or comparisons to forward models have been used 
to understand broad trends in the brown dwarf popula¬ 
tion. Such analyses have led to the L-T-Y classification 
schemes as well as widely accepted spectral and photo¬ 
metric indicators of gravity and metallicity (Kirkpatrick 
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2005; Gushing et al. 2011). However such approaches 
have important limitations. The physics and chem¬ 
istry of these dense molecular atmospheres is complex 
and subtle processes can significantly influcence thermal 
spectra. 

The traditional method for interpreting brown dwarf 
spectra is to construct sets of sophisticated cou¬ 
pled radiative-convective-chemical equilibrium models to 
which the data can be compared (Marley et al. 1996; Al¬ 
lard et al. 1996; Tsuji et al. 1996; Burrows et al. 2001). 
Systematic comparisons of models to data constrain at¬ 
mospheric composition, temperature, and cloud proper¬ 
ties (e.g., Saumon et al. 2006; Gushing et al. 2008; Rice 
et al. 2010). While these grid modeling studies have un¬ 
questionably advanced our understanding of brown dwarf 
atmospheres, there is still much that we don’t under¬ 
stand. For instance, can dynamical processes in brown 
dwarf atmospheres cause deviations from radiative con¬ 
vective equilibrium? How do their molecular composi- 
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tions vary with altitude? Can we measure the molecular 
abundances rather than assume them based upon chem¬ 
ical equilibrium or disequilibrium? Can their elemental 
abundances deviate from the often assumed solar com¬ 
position? It is difficult to answer such questions from 
comparison solely with forward models, as not all pro¬ 
cesses can be modeled with fidelity and, in any case, it 
is not always obvious which processes are responsible for 
given deviations of models from data. 

Line et al. (2014b) presented a novel approach for an¬ 
swering the above questions. Rather then rely on grid 
models that only constrain a few basic parameters, they 
used atmospheric retrieval methods common in Earth 
and planetary atmosphere studies to invert brown dwarf 
spectra for the temperature structures and abundances. 
Using such approaches allows the maximum amount of 
information to be extracted from brown dwarf spec¬ 
tra. Line et al. (2014b) found for, G1570D, with few 
assumptions about the nature of the chemistry or the 
temperature-pressure profile, that the retrieved quanti¬ 
ties were consistent with previous grid modeling studies 
(Saumon et al. 2006), though with some minor devia¬ 
tions such as a more isothermal upper atmosphere and 
a depleted ammonia abundance. While individual ob¬ 
jects themselves are interesting, the key to understand¬ 
ing the underlying physical processes in brown dwarf at¬ 
mospheres (or any class of atmospheres) requires an in¬ 
vestigation of a large population of objects. With large 
populations one can identify trends and correlations of 
various physical parameters that can lead to insight into 
new phenomena. 

We first aim to improve upon the techniques presented 
in Line et al. (2014a) and to validate our improved ap¬ 
proach against benchmark brown dwarfs. Then in a fol¬ 
lowup paper we will apply our new approach to a small 
sample of T-dwarfs. In that sample we will identify 
trends within our retrieved results such as metallicity vs. 
gravity, molecular abundances vs. effective temperature 
etc., and also how empirical spectral properties such as 
color or other indices correlate with the retrieved physi¬ 
cal parameters. Preliminarily, we choose mid- to late-T 
dwarfs as they, based upon previous investigations, ap¬ 
pear to be largely free of thick silicate clouds that can 
complicate the interpretation of L and early-T spectra 
(Kirpatrick 2005). 

In this manuscript we present the upgraded retrieval 
approach and apply it on two benchmark brown dwarfs, 
HD3651B and G1570D, using the the SpeX Prism Library 
(Burgasser et al. 2006a). In ® we present a modified 
Bayesian retrieval approach anoa novel approach for in¬ 
verting for temperature structures. In ^we present our 
retrieved results and how they are impacted by various 
assumptions. Finally, in Qwe compare our retrieved val¬ 
ues to benchmark properties such as the age and metal- 
licity. Additionally, we present a detailed stellar abun¬ 
dance analysis of G1570A and HD3651A and derive their 
carbon-to-oxygen ratios in order to add another dimen¬ 
sion to the benchmark comparison. 

2. METHODS 

Here we give a brief review of current state of knowl¬ 
edge of G1570D and HD651B, describe the data we use, 
the forward radiative transfer model, and the inverse 
methods to retrieve temperature and abundance infor¬ 


mation. Building upon Line et al. (2014b), we have 
made significant upgrades to our methodology in terms 
of the forward model, Bayesian estimator, and treatment 
of the data. We highlight those differences where appli¬ 
cable. 

2.1. Current State of Knowledge on G1570D & 
HD3651B 

The targets of this study have been selected for their 
suitability as robust test cases for validating our ap¬ 
proach. Both objects are wide-orbit common proper mo¬ 
tion companions to stars, allowing us to check our derived 
properties for consistency with those found for their stel¬ 
lar hosts by other routes. For this reason, such systems 
are often referred to as benchmarks, although the qual¬ 
ity and context of the available constraints varies widely 
depending the mass and evolutionary phase of the stellar 
primary (e.g., Pinfield et al. 2006). 

G1570D was the first T-dwarf companion to a star iden¬ 
tified following the prototypical T-dwarf G1229B (Bur¬ 
gasser et al. 2006a). It is a wide component in a hi¬ 
erarchical quadruple system, whose inner components 
are an M1V+M3V spectroscopic binary (G1570B and 
G) and a K4V primary (G1570A) (Gleise 1969; Duquen- 
nory & Mayor 1988; Mariotti et al. 1990; Foreville et 
al. 1999), from which G1570D lies at a projected sepa¬ 
ration of 1525 ± 25 AU (Burgasser et al. 2000). G1570D 
has been subject to a number of grid-model fitting stud¬ 
ies, which have to varying extents used the primary star 
to restrict the parameter space available for the models 
(e.g., Geballe et al. 2001; Saumon et al. 2006; Legett et 
al. 2007; Saumon et al. 2012) and has been used as an 
anchor point for applying trends seen in self-consistent 
grid models to estimate parameters of the wider T dwarf 
population (Burgasser et al. 2006b). 

Liu, Leggett & Ghiu (2007) used a variety of stellar 
age indicators for G1570A to constrain the age of the 
system to the range 1-5 Gyr, whilst Saumon et a. (2006) 
collated literature values to estimate its metallicity as 
[Fe/H] = 0.09 ± 0.04. In the Appendix we present a new 
measurement of [Fe/H] = —0.05 ± 0.17. 

HD3651B was identified as a wide common proper mo¬ 
tion companion to the planet-hosting KOV star HD3651A 
with a projected separation of 480 AU (Murgrauer et al. 
2006; Luhman et al. 2007; Liu, Leggett & Chiu 2007). 
Like G1570D, it has been the subject of a number of 
spectroscopic studies that have been constrained by the 
properties of the primary star (e.g., Leggett et al. 2007; 
Burgasser 2007). Liu, Leggett & Chiu (2007) reviewed 
X-ray luminosity, chromospheric and rotation-based age 
indicators for HD3651A, and found the target lies in the 
unreliable “older” tail of each of these diagnostics. Like 
them, we adopt the isochronal age range of 3-12Gyrs 
from Velenti & Fischer (2005). As an exoplanet host star, 
HD3651A has been the subject of several recent compo¬ 
sition studies. The determinations of [Fe/H] are consis¬ 
tently super-Solar, ranging from [Fe/H] = +0.12 ± 0.04 
(Santos et al. 2004) to [Fe/H] = +0.19 ± 0.03 (Ghezzi 
et al. 2010). In the Appendix we have further devel¬ 
oped these targets potential to contribute to our under¬ 
standing its substellar companions by compiling detailed 
abundance measurements from the literature, and mea¬ 
suring new abundances for both. Most significantly for 
this study we present new determinations of the C to O 
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ratios for both primary stars. 

2 . 2 . Data 

We use the data within the SpeX Prism Library (Bur- 
gasser et al. 2006) to perform our analysis on our tar¬ 
get objects. Since a given SpeX spectrum is continuous, 
we avoid having to consider various instrumental sys- 
tematics (e.g., as had to be done in Line et al. 2014b) 
and subsequent impact on their interpretation that come 
along with having to stitch spectra from multiple instru¬ 
ments together. We use the SpeX data taken in the SXD 
mode which cover wavelengths between 0.8 - 2.5 /am with 
a wavelength dependent resolving power ranging from 
87 to 300. The spectra within the SpeX SXD library 
are oversampled relative to a spectral resolution element, 
and are therefore each pixel is not an independent sam¬ 
ple. Using the full oversampled data would result in over 
constrained results. We therefore sample every few pix¬ 
els. We choose the sampling length based upon the auto- 
corrletion length scale of the residuals of a typical model 
fit to the data. This is 2.7 pixels. We thus take every 
other 3rd pixel to be statistically independent. An alter¬ 
native approach would be to model the covariant error 
structure of the oversampled data through a Guassian 
process (e.g., Czlecka et al. 2014). 

The native format of the spectral fluxes and error bars 
within the library are normalized spectra. In order to 
perform the subsequent analysis, these spectra and er¬ 
ror bars must be converted into physical units via pho¬ 
tometric calibration. The SpeX database provides the 
2MASS photometric J, H, and K-band magnitudes for 
each object. We convert the 2 MASS magnitudes into 
MKS flux units (W m^ m“^) using the Spitzer Science 
Center Magnitude/Flux density converter R which uses 
the zero point fluxes described in Cohen, vVheaten & 
Megeath (2003) . The H-band flux is used to derive 
the final flux-calibrated spectrum just as in Saumon et 
al. (2006) and Liu, Leggett & Chiu (2007). Figure 
shows the flux calibrated spectra. Unfortunately, the 
calibrated spectra are different by some scale factor that 
is larger than the quoted photometric uncertainty, de¬ 
pending upon which photometric point is used to cal¬ 
ibrate. Our retrieval model (discussed below) includes 
a scaling factor as a free parameter so these differences 
are not critical to our analysis (with the exception of the 
derived spectroscopic radius, see the Results section for 
more on this). We do note, however, that better pho¬ 
tometry, namely (MKO) exists for these objects. Fur¬ 
thermore, Stephens & Leggett (2004) suggest that the 
2 MASS photometry is not the most accurate due to the 
shape of the 2 MASS filter profiles with respect to telluric 
transmittance. They also provide correction factors for 
the 2 MASS photometry based upon the more accurate 
MKO photometry. However, we choose to use the un¬ 
corrected 2 MASS photometry as they shall provide the 
most conservative impact of inconsistencies in photome¬ 
try on the derived quantities. 

2.3. Forward Radiative Transfer Model 

The forward radiative model is a derivative of the 
CHIMERA forward model (Line et al. 2013a; 20I4a,b) 

^ http://ssc.spitzer.caltech.edu/warmmission 
/ propkit / pet / magto jy 
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Fig. 1. — Spectral Calibration Process on the two objects. The 
photometric (circles) J (blue), H (green), and K (red) points are 
used to calibrate the normalized SpeX spectra. The spectra are 
calibrated by integrating over the corresponding filter profile and 
rescaling the normalized spectrum to match the photometric fluxes. 
Each photometric point gives a different calibrated spectrum-blue 
for J-band, green for H-band, and red for K-band. The H-band 
calibrated spectrum with the error bars is shown with gray dia¬ 
monds and is what we use to perform our analysis. Note that the 
differences between each band integrated spectrum greatly exceeds 
that of the photometric uncertainty. 

which computes the upwelling I-dimensional disk in¬ 
tegrated thermal emission spectrum given the molecu¬ 
lar abundances, temperature-pressure (TP) profile, and 
gravity g. Near infrared spectra of late-T’s are typi¬ 
cally dominated by strong absorption features from wa¬ 
ter, methane, and alkali metals and little if any obvi¬ 
ous absorption exists due to other gases. We include 
constant-with-altitude volume (molar) mixing ratios for 
H2O, CH 4 , CO, CO2, NH 3 , H2S, and alkali opacities. 
These are the species known to be found in cool dwarf 
atmospheres that have spectral signatures in the near 
infrared. The alkali opacities include only sodium and 
potassium and are treated as only one free parameter 
with their ratio assumed to be solar. Hydrogen/helium 
in solar ratio is assumed to make up the remainder of the 
gas. 
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Th e TP profile is also included as a free parameter 
(see ^2.4.2). The TP profile is partitioned into 15 evenly 
spaced slabs (or knots) in log pressure between 315 bar 
and 1 mbar. Because a 15 layer atmosphere does not have 
enough vertical resolution to accurately compute fluxes, 
the 15 level profile is spline interpolated to a finer 70 level 
grid before before being passed to the radiative transfer. 
Using fewer TP points permits swifter convergence of the 
Bayesian estimator. 

In addition to the TP profile and gas mixing ratios, 
we include gravity as a free parameter. Gravity con¬ 
trols the column optical depth. Most of the opacity 
database is drawn from Freedman et al. (2014) and ref¬ 
erences therein as in Line et al. (2014b), but we have 
also incorporated the most up-to-date methane line list 
(Yurchenko et al. 2014) using the line broadening coeffi¬ 
cients from Margolis (1996). Instead of using line-by-line 
or correlated-K we simply sample the hi-resolution cross 
sections at 1 cm“^ resolution (see Sharp & Burrows 2007 
section 2) which is more than sufficient for moderate res¬ 
olution spectra. 

Finally, the high resolution spectra are convolved with 
a wavelength dependent Gaussian instrumental profile 
that reflects the wavelength dependent resolving power, 
and then binned to the data wavelength grid in order for 
direct data-model comparison. 

Figure shows the sensitivity of a model spectrum 
at SpeX resolutions to the various parameters. Many 
of these parameters are sensitive to similar wavelengths. 
This will result in correlations/degeneracies amongst the 
gases and the temperature at different levels in the at¬ 
mosphere. 

All of the aforementioned parameters (23 of them) con¬ 
trol the flux at the top of a brown dwarf atmosphere. We 
also include additional “systematic” parameters that fa¬ 
cilitate the direct comparison of the model to the data. 
These parameters account for the radius to distance ratio 
(and implicitly the flux calibration), uncertainties in the 
wavelength calibration, underestimation of the spectral 
error bars, and a smoothing parameter for the temper¬ 
ature profile-for a total of 27 free parameters. These 
“systematic” parameters will be discussed in more detail 
in the following section. A list of all of the parameters 
and a brief description of each is presented in Table 
The model has many parameters, but this allows us to 
make very few implicit assumptions about the nature of 
the atmosphere. Unconstrained parameters will simply 
appear unconstrained. The beauty of modern retrieval 
approaches (below) is that they can accommodate nu¬ 
merous parameters and fully account for all of the cor¬ 
relations amongst them. A larger number of parameters 
will of course result in more conservative estimate of the 
uncertainties through marginalization. 


TABLE 1 

Parameters in the forward model. 


Parameter 

Description 

log/z 

log of the volume mixing ratios of H 2 O, 

CH 4 , CO, CO 2 , NH 3 , H 2 S, and alkali (Na±K) 

logs' 

log gravity [cms“^] 

(R/D)2 

radius-to-distance scaling (Rj/pc) 

T,- 

temperature at 15 pressure levels (K) 

AA 

uncertainty in wavelength calibration (nm) 

b 

error bar inflation exponent (equation 

7 

TP profile smoothness hyperparameter 



1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 

A (/^m) 


2.4. Retrieval Model 
2.4.1. Bayesian Implementation 

The retrieval model is the Bayesian engine that op¬ 
timizes the forward model to fit the data. We use the 
Markov chain Monte Garlo approach implemented with 
afhne-invariant ensemble sampler, EM GEE (Eoreman- 
Mackey et al. 2013). This is a significant advancement 
over the optimal estimation and bootstrap Monte Garlo 
approaches used in Line et al. (2014b) as we are now able 


Fig. 2.— Sensitivity of the spectrum to various parameters. This 
is a synthetic spectrum with an effective temperature (equivalent 
black body flux integrated over the whole spectrum) of 700 K and 
a log^f of 5 and with purely thermochemical equilibrium composi¬ 
tion. The red regions represent a change in the spectrum due to a 
perturbation of each of the parameters. For H 2 O, CH 4 , NH 3 , H 2 S, 
alkali, this perturbation is ± 0.5 dex (where 1 dex is one increment 
in log space, or 1 order of magnitude) in number mixing ratio from 
the thermochemical abundance value. For CO the perturbation is 
+2 dex, CO 2 , +6 dex, log^fT 0.1 dex, and the temperatures are 
perturbed at each level by ±50 K. 
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to make fewer a priori assumptions about the smoothness 
of the temperature profile or the Gaussian shape of the 
parameter uncertainties. EMCEE requires only a func¬ 
tional form for the log of the posterior probability to 
perform the optimization. The posterior probability is a 
combination of the likelihood and the prior described as 
follows. Starting from Bayes theorem 


p(x|y) 


£(y|x)p(x) 

E 


( 1 ) 


where x is the parameter vector described in §2.3| and 
y is the data vector-in our case the spectrum, p(x|y) 
is the posterior probability distribution, >C(y|x) is the 
likelihood distribution which penalizes poor fits to the 
data, p(x) is the prior which represents any external con¬ 
straints, and E is a normalization factor known as the 
evidence, or marginal likelihood, which is required for 
Bayesian model comparison but not for parameter esti¬ 
mation. We use the following log-likelihood function: 


ln«yW = -iX:^ 


i=l 


TW)!_lln(27rs?) (2) 


Here, the index i denotes the data point, in our case 
some property at a single wavelength bin, y is the mea¬ 
sured flux, E(x) is t he m odeled flux that comes out of 
the forward model ( ^2.3[ ), and s is the data error given 
by 

sf = af + lo" (3) 

where a is the measured error for the data point 

and 6 is a free parameter. Differing from Line et al. 
(2014b), we modify the standard error on the data point 
by the factor 10 ^ to account for underestimated uncer¬ 
tainties and/or unknown missing forward model physics 
(Eoreman-Mackey et al. 2013, Hogg et al. 2010, Tremain 
et al. 2002), e.g., imperfect fits. This results in a more 
generous estimate of the parameter uncertainties. Note 
that this is similar to inflating the error bars post-facto in 
order to achieve reduced chi-squares of unity, except that 
this approach is more formal because uncertainties in 
this parameter are properly marginalized into the other 
relevant parameters. Generally, the factor 10 ^ takes on 
values that fall between the minimum and maximum of 
the square of the data uncertainties. The first term inside 
the summation in equationis the familiar “chi-square”. 
This term penalizes large residuals. The second term 
in the summation is the Gaussian normalization factor 
that is normally excluded from standard fitting routines 
due to the unchanging data errors. Because the data er¬ 
rors include the free parameter b this normalization can 
change, and hence has to be taken into account. Re¬ 
ally, the purpose of this term is to provide a balance 
for the error bar inflation parameter to prevent it from 
approaching infinity. 

The prior, p(x), can be broken up into several pieces 
as 

p(x|y) = p(T)p(x')p(7) (4) 

where p(x') is the prior on the log of the gas mixing 
ratios, the instrumental parameters, gravity, and the 
radius-to-distance scaling while p(T) and ^(y) are the 
temperature profile priors. The parameter 7 is the rela¬ 


tive weighting of the temperature prior (see equation]^ 

. The prior details are shown in Table 

Because we have both measured fluxes and parallaxes 
for each object, we are able to calculate the “photomet¬ 
ric” radius. If we can measure both radius and a gravity 
we can then constrain the mass. We know for brown 
dwarfs the mass cannot exceed ^ 75 — 80Mj (e.g.. Bur¬ 
rows et al. 2001). Therefore, rather than place individual 
priors on the radius and gravity we enforce a prior on the 
derived mass to fall between the physical plausible values 
of 1 and 80 Mj. This constraint prevents the retrieved 
radii and gravities from entering an unphysical region of 
parameter space. 

2.4.2. A Novel Temperature-Pressure Profile Retrieval 
Approaeh 

We present our novel method for retrieving tempera¬ 
ture profiles in atmospheres. A common issue in plan¬ 
etary atmospheric retrievals is how to parameterize the 
temperature profile in an atmosphere. There are two 
philosophies. One philosophy, mostly used in the exo¬ 
planet atmosphere community when the data is sparse, 
is to parameterize the atmosphere with some analytic 
function that can be described by a small set of parame¬ 
ters (e.g., Madhusudhan & Seager 2009; Line et al. 2012 ; 
2013 Benneke & Seager 2012 ). This is advantageous 
as the entire TP structure can be controlled by just a 
few simple parameters. It is disadvantageous because it 
is relying on a parameterization to infer the tempera¬ 
ture structure. This could potentially result in biases in 
the retrieved profiles. Eor instance, if one were to use 
the simple Eddington approximation for the tempera¬ 
ture profile, there would be one free parameter-the mean 
opacity. While just one free parameter is ideal, we know 
that the Eddington approximation is a poor approxima¬ 
tion for brown dwarf atmospheres because the true opac¬ 
ity is not constant with pressure, nor is it gray. While 
parameterizations are appealing in their simplicity, they 
can often times be too much of an oversimplification of 
the physics, and are thus not appropriate. 

The classic planetary science approach, the other ex¬ 
treme, is to retrieve the temperature at each model 
layer in the atmosphere (e.g., Rodgers 2000; Irwin et 
al. 2008; Lee et al. 2012 ) within an optimal estimation 
framework. Eor typical model atmosphere grids this can 
be anywhere between 50-100 independent temperature- 
pressure points. This was the approach used in Line et 
al. (2014b). Because many of the atmospheric levels are 
degenerate, the retrieved profiles often result in unphys¬ 
ical oscillations, or ringing (Rogers 2000). The standard 
remedy to this problem is to implement some a priori co- 
variance matrix, or a smoothing kernel (or Tikonov Reg¬ 
ularization), given some set smoothing length scale and 
a priori width (Irwin et al. 2008). While this reduces 
wild oscillations, this smoothing length scale and width 
must be chosen a priori and cannot change during the 
course of a retrieval. Eurthermore, these values are often 
case specific and must be tuned by hand. Our novel ap¬ 
proach remedies all of the aforementioned issues, and can 
be readily implemented within a Bayesian framework. 

We borrow much of this work from the non-parametric 
regression literature (Lang & Brezger 2004; Rahman 
2005; Jullion & Lambert 2007). The goal is to allow 
flexibility to fit for each of the independent temperature- 
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pressure points while preserving smoothness, without 
having to a priori set the degree of smoothness. The 
best way to do that is by penalizing the second derivative 
of the temperature structure. The second derivative is 
the “roughness” of a function. Temperature vectors with 
wild, unphysical oscillations, will have large summed sec¬ 
ond derivatives. These types of roughness-penalized non- 
parametric polynomial fits are known as P-splines. The 
degree to which this roughness is penalized is included 
as a free parameter ( 7 ). This variable smoothing is im¬ 
plemented as, 

1 ^ 1 
Inp(T) = - 2^* +Ti_i)2 - - ln(27r7) (5) 

Inside the sum is the discrete second derivative of the 
temperature profile at each level, i weighted by 7 . Based 
on experimentation (Lang & Brezger 2004; Rahman 
2005; Jullion & Lambert 2007) the hyperprior on 7 
should take the form of an inverse gamma distribution 
with the properties shown in Table A variable 7 al¬ 
lows the data to dictate the degree m smoothing. If the 
data warrants little smoothing and there truly are os¬ 
cillations in the TP profile, then 7 will be large, lending 
little weight to the smoothing prior. When the data does 
not justify rough TP profiles, 7 will be small resulti ng in 
a larger penalty to rough profiles. As mentioned in §2.3[ 
15-knots, or anchor points are used in our TP prohle. 
This number is somewhat arbitrary as the number and 
location of the knots should not matter within this frame¬ 
work just so long as there are enough evenly spaced knots 
to sufficiently resolve potential structure (Eilers & Marx 
1996). We have tested higher numbers of knots (30 vs. 
15) and have indeed found little sensitivity to the choice. 

The combined log-likelihood and log-prior are readily 
implemented within the EMCEE sampler. The MCMC 
is initialized with a tight Gaussian ball with 8 chains 
or walkers (Eoreman-Mackey et al. 2013) per parameter 
about an educated starting guess in order to minimize 
burn in time (the time it takes for the MCMC sampler 
to locate the posterior distribution). We note that the 
final results are insensitive to the initial starting point; 
poor starting points result in longer burn in times. Con¬ 
vergence is monitored with the Gelman-Rubin statistic. 
The statistics are summarized by drawing thousands of 
random points from the last ^ 5 — 10% (which we take 
to be the posterior) of the ensemble of chains which gen¬ 
erally entail 500-1000 independent samples (as dictated 
by the autocorrelation length scale). In the following 
sections we describe the physical properties of the two 
benchmark late T-dwarfs evaluated within this frame¬ 
work. 

3. RETRIEVAL RESULTS 

The fiducial retrieval results are summarized in Eig- 
ures[^and[^ The top row of Eigure shows an ensem¬ 
ble 01 thousands of fits derived from the posterior and 
their residuals. We also summarize the retrieved log^f 
and Teff. The residuals appear to be minimal and ran¬ 
dom, with exception of the peaks near Y and J band. 
The slight misfit at these wavelengths could potentially 
be due to uncertainties in the alkali (sodium, potassium) 
cross sections, which have the least certain cross-sections 
of the opacities that absorb at these wavelengths. We 


TABLE 2 

Summary of priors for each of the parameters. 


Parameter 

Prior 

log A 

uniform-in-lo^with log A ^ “12, E^A ^ 1 

(R/D)2, log^ 

uniform, constrained by 1 <qI^ jG < 80Mj 

Ti 

see equation 

AA 

uniform (-10 - 10 nm) 

b 

uniform, 0.01 x < 10^ < 100 x max{cr‘^) 

7 

Inverse Gamma (r( 7 ; a,/3)), a = l,d = 5 X 10“^ 


^We note that uniform-in-log priors may cause issues for a larger 
parameter set. For larger parameter sets a Dirichlet prior should 
be used 


note that these residuals are much smaller than can be 
obtained by typical grid model fits suggesting that the 
additional parameters we include in our model are re¬ 
quired to produce these better fits. 

The retrieved temperature profiles (bottom row, Eig- 
ure 1 ^ for each object are summarized with a median, 
1 , and 2(7 credibility region from the ensemble of thou¬ 
sands of randomly drawn TP profiles from the poste¬ 
rior. The retrieved TP profiles appear to be consistent 
with physically based expectations. In each panel we 
show for comparison a self-consistent grid model profile 
(black, from the grid models of Saumon & Marley 2008) 
corresponding to the median retrieved effective tempera¬ 
ture and gravity. The self-consistent grid models assume 
cloud-free 1 -dimensional radiative convective equilibrium 
with a solar composition atmosphere in thermochemical 
equilibrium. The agreement is astounding and this is a 
point we would like to stress. The self-consistent grid 
model profiles generally fall well within the 2(7 credibil¬ 
ity region. We have made no assumptions about the 
nature of the TP profiles other than smoothness, the 
degree of which was allowed to vary. We also tested dif¬ 
ferent starting guesses (e.g., isothermal) and the results 
are no different. This suggests that these are the actual, 
true, temperature-pressure profiles in these atmospheres 
and that the assumption of one-dimensional radiative- 
convection is sufficient. However, there is some diver¬ 
gence at pressure levels less than about 1 bar and greater 
than a few 10s of bars. The retrieved profiles tend to be¬ 
come more isothermal near the top of the atmosphere 
than the radiative convective models, though this diver¬ 
gence is much smaller than the width of the confidence 
intervals and is therefore not significant (as constrained 
by SpeX data alone). We also note that a recent paper 
by Tremblin et al. (2015) predicts that condensation- 
induced fingering convection can result in a cooler deep 
atmosphere than expected from standard dry convection 
assumptions, consistent with what we are finding but not 
conclusive. Higher resolution data with more vertical 
resolution and altitude range or longer wavelength data 
will (and have, as in Line et al. 2014b) provide better 
constraints to the TP profile that will allow us to further 
test deviations from the standard radiative-convective as¬ 
sumptions. In this investigation, we purposefully avoid 
combining different datasets in this investigation due to 
the introduction of additional systematics that come with 
multiple data sets. 

Eigurej^ summarizes the posterior for the physical at¬ 
mospheric parameters. The effective temperature, ra¬ 
dius, mass, metallicity ([Ee/H]), and C/0 ratio are all 
derived quantities and were not directly retrieved. The 
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Fig. 3. — Spectra (top row) and retrieved temperature profiles (bottom row). For the two objects we show the H-band calibrated SpeX 
data as the diamonds with error bars, a summary of thousands of model spectra generated from the posterior and their residuals (median 
in blue, Icr spread in red), and their spectral type and bulk properties. The bottom row summarizes thousands of temperature profiles 
drawn from the posteriors for each object (median in blue, Icr spread in red, 2cr spread in pink). The black temperature profile shown for 
each object is a representative self-consistent grid model (Marley et al. 2008) interpolated to the quoted log^f and Tgff to demonstrate that 
our retrieved profiles are physical and are consistent with 1-D radiative convective equilibrium. 


effective temperature distribution is obtained by comput¬ 
ing the bolometric flux ( 1-20 /im) over thousands of spec¬ 
tra generated from the posterior. The radius is derived 
from the retrieved spectral scaling factor, (R/D)^, given 
the distance. The distance and photometric uncertain¬ 
ties are formally propagated into the radius uncertainty 
via Monte Carlo error propagation. The mass is derived 
from the retrieved gravity and radius. We also reiterate 
that the radius, and hence the derived mass, are very 
sensitive to the accuracy (e.g., missing systematics not 
accounted for in the quoted photometric errors) of the 
photometry. We discuss the impact of the photom etric 
calibration on the retrieved quantities further in §3.3[ 
The metallicity is derived by summing up the molecular 
mixing ratios for each species weighted by the number of 
metal atoms divided by the abundance of hydrogen and 
then comparing that to the sum of solar metals relative 
to hydrogen. The C to O ratio is computed by divid¬ 
ing the sum of the carbon bearing species by the oxygen 
bearing species appropriately weighted by the number of 


carbon/oxygen atoms in each species. 

The marginalized probabilities for each parameter (the 
histograms) are shown along the diagonals for each ob¬ 
ject. For both objects the H 2 O, CH 4 , NH 3 and Na+K 
are constrained ( 68 % confidence) to better than 0.3 dex 
(or a factor of 2). For comparison the best constraints 
we have on gases in exoplanet atmospheres are to within 
a factor of ^ 10 (e.g., Kreidberg et al. 2014). 

Perhaps the most surprising finding is the tight con¬ 
straint on ammonia for each object. The robustness of 
the am monia constraints in the near-IR are discussed 
in §3.1[ Additionally there are strong correlations of 
log^f with spectrally prominent absorbers, H 2 O, CH 4 , 
and NH 3 . As expected for any atmosphere in hydro¬ 
static equilibrium, there is a positive correlation between 
metallicity and gravity (this can also be seen in the 
metallicity vs. gravity panels in Figure [^. As grav¬ 
ity increases, the optical depth at a given layer in the 
atmosphere decreases (r = where n is the opac¬ 

ity), so the opacity must increase to maintain that same 
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Fig. 4.— Summary of the posterior for the relevant parameters. The stair-step plot on the top right is for G1570D, and the bottom left 
for HD3651B. These show the marginalized posteriors (ID histograms) along the diagonals and the parameter correlations (2D histograms). 
The parameters for each object are on the same scale and so can be directly compared. The dashed lines in the ID histograms are the 
16, 50, and 84 percentiles. The width between the 16 and 84 percentiles represents the 68% confidence interval. For each parameter, the 
median and plus/minus Icr values are shown just above (HD3651B) and below (G1570D) the histograms. 
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optical depth. The strong absorbers, H 2 O, CH 4 , and 
NH 3 , are also correlated with each other and with tem¬ 
perature (not shown). These correlations are positive, 
which seems backwards for overlapping absorbers, but 
it is because as one absorber increases, the temperature 
first increases to maintain that same flux at the wave¬ 
length of that absorber, resulting in a higher flux at a 
different wavelength where another absorber is present. 
This absorber must increase in abundance to suppress 
the flux at this different wavelength. The other gases, 
CO, CO 2 , and H 2 S, are largely unconstrained; only up¬ 
per limits can be obtained. This is mainly because of 
their relatively low thermochemical abundances (despite 
potential vertical mixing) and relatively weak bands in 
the near infrared. We would expect CO and CH 4 to flip 
roles for objects hotter than ^1100 K. 

In order to check whether or not the retrieved abun¬ 
dances are chemically realistic, we compare them to ther¬ 
mochemical equilibrium models. To do this, we use the 
NASA Chemical Equilibrium with Applications model 
(CEA, Gordon & McBride 1996) recently used for exo¬ 
planet atmosphere studies (Line et al. 2010; Moses et al. 
2011; Line et al. 2011). CEA only requires the local tem¬ 
perature, pressure, and elemental abundances at given 
model layer in order to compute the equilibrium abun¬ 
dances. Equilibrium condensation chemistry (no rain- 
out) is included but the code has a difficult time with 
enstatite (MgSiOs) condensation. Eor this reason we do 
not include magnesium or silicon species, but account for 
the depletion of oxygen due to presumed cloud formation 
in the deep atmosphere by removing 3.28 oxygen atoms 
for every silicon atom (Burrows & Sharp 1999). We have 
not accounted for perturbations to equilibrium chemistry 
such as horizontal or vertical mixing. 

Eigure compares the retrieved results for the well 
constraint species to thermochemical equilibrium abun¬ 
dances along the median temperature profile. In princi¬ 
ple there will be a spread, albeit minor, due to uncer¬ 
tainties in the temperature profile, but since the goal is 
to just check for realism in the abundances we ignore 
such a spread. Eor each object we show two cases of 
equilibrium abundances. The first is for solar elemental 
composition (solid lines) and the second is a by-hand fit 
of the intrinsic metallicity and C to O ratio to the re¬ 
trieved quantities. We stress that the intrinsic metallic¬ 
ity and C/0 are different from the atmospheric metallic¬ 
ity and C/0. Condensate processes can deplete oxygen or 
other species in the atmosphere resulting in atmospheric 
metallicity and C/O’s that differ from the intrinsic or 
bulk values. 

We find for both G1570D and HD3651B that the as¬ 
sumption of intrinsic solar elemental abundances over¬ 
estimates the retrieved water and methane abundances 
but appears to do a good job for the alkali and ammo¬ 
nia mixing ratios. However, by hand-tuning the bulk 
metallicity and C/0 in the thermochemical model we 
can better match the retrieved water and methane abun¬ 
dances. This is not a rigorous “fit” to the chemistry by 
any means, but simply an attempt to show that we have 
retrieved chemically plausible molecular abundances. A 
perhaps more rigorous approach for obtaining the al¬ 
lowed ranges of C to O ratios and metallicities would 
be to perform a “retrieval on the retrieval” where by the 
chemical model would fit the retrieved molecular abun¬ 


dances within a Bayesian framework. This is currently 
beyond the scope if this work. We also note that the 
ammonia thermochemical profiles agree well with our 
column averaged uniform-with-altitude retrieved values. 
Saumon et al. (2006) suggest that quenching of ammo¬ 
nia due to vertical mixing occurs in the deep atmosphere 
near a temperature of ^2200 K. Such temperatures occur 
at the deepest pressures (several hundred bars) on our 
retrieved profiles. If ammonia indeed quenches at these 
deep levels then the ammonia profile would be nearly 
constant with altitude well within our retrieved range 
for both objects (Eigure [^. 

Another surprising find is that the retrieved wa¬ 
ter/methane abundance is lower (water/methane ^ 0 . 8 - 
0.9) than what one may expect for a solar composition 
atmosphere. Typically water is more abundant than 
methane by a factor of ^1.5 at solar composition (see 
Saumon et al. 2006). This suggests a super-solar (greater 
than 0.5) atmospheric carbon-to-oxygen ratio. Account¬ 
ing for the draw down of O due to silicate condensation, 
the methane/water abundance shown in Saumon et al. 
(2006) for G1570D suggests an atmospheric C/0 of 0.63. 
Our retrieved atmospheric carbon-to-oxygen ratios for 
both objects are higher than one. The inferred, by-hand 
intrinsic C to O ratios are less than unity but are still 
slightly higher than solar (see ^for a comparison to the 
host star values). We note that, in our previous study. 
Line et al. (2014b), we found a fairly low (^ 0 . 2 ) carbon- 
to-oxygen ratio for G1570D. These differences are likely 
due to the inclusion of the Akari and Spitzer Infrared 
Spectrometer data and the treatment of the systemat- 
ics between them. Our current investigation is more 
straightforward as we only focus on the SpeX data set, 
and thus do not have to worry about potential biases due 
to differing unaccounted for systematics amongst differ¬ 
ent datasets. In Qwe show that our thermochemically 
self-consistent intrinsic carbon-to-oxygen ratios for the 
brown dwarf’s are consistent with the those derived from 
the stellar primaries. 

3.1. Ammonia in the Near-IR 

One of the more remarkable findings in this investi¬ 
gation is the strong evidence for the presence of ammo¬ 
nia in low resolution near infrared speetra. Ammonia at 
longer wavelengths in T-dwarfs is not new. Cushing et 
al. (2008) and Saumon et al. (2006) convincingly demon¬ 
strated the presence of ammonia in G1570D using the 
strong 9.6 pm band in the Spitzer IRS data. However, 
spectroscopic features of ammonia are not expected to 
present themselves below 2.5 microns until the Y-dwarfs 
(Kirpatrick et al. 2005) unless observed at high resolu¬ 
tion (Canty et al. 2015). Therefore, we were surprised to 
find how strongly ammonia can be constrained (e.g., an 
actual bounded limit as opposed to an upper limit only) 
with /ore-resolution near infrared data in both objects 
despite the lack of obvious spectral features in this wave¬ 
length range. How are we to believe that our constraint 
is real? We show three lines of evidence supporting our 
strong ammonia constraint. 

One line of evidence comes from a standard Bayesian 
hypothesis testing procedure. Such procedures deter¬ 
mine whether or not a parameter within nested mod¬ 
els is justified given the data (e.g., Trotta et al. 2008). 
A commonly used approach is to compare the Bayesian 
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Fig. 5. — Comparison of the retrieved values (shaded boxes) of the well constrained molecules with their expected thermochemical 
equilibrium abundances along the median temperature profile. The solid curves are the thermocemical equilibrium abundances for solar 
composition while the dashed curves are the thermochemical equilibrium abundances for the specified C/O and metallicity. This shows 
that the retrieved abundances are thermochemically consistent. 


information criterion (BIC) between a model with and 
without a particular parameter. Parameters that pro¬ 
vide better fits to the data and produce a delta-chi-square 
that is greater than the penalty of adding that additional 
parameter are justified. However, the BIC is a truncated 
Laplace approximation to the full Bayesian evidence, or 
marginal likelihood (Kass & Raferty 1995). We do not 
use the BIC here, rather we compute the full Bayesian ev¬ 
idence. This is done by numerically integrating over the 
entire posterior using the approach described in Wein¬ 
berg et al. ( 2012 ) (see also Swain, Line & Deroo 2014 for 
an application to exoplanet spectra). By computing the 
evidence of the full model that contains all parameters 
to the one that removes ammonia, we can obtain a Bayes 
factor. Bayes factors greater than one suggest that the 
model containing the parameter in question is favored, 
while Bayes factors less than one suggest otherwise. A 
Bayes factor can then be converted into a confidence of 
detection (Trotta 2008). Table shows the Bayes factors 
and the corresponding detection significances for three 
different nested models each removing only one gas (NH 3 , 
H 2 S, or H 2 O) from the full model in Tablef^ 

We show the detection significances for H 2 O as an ex¬ 
ample of a gas that is visibly obvious in the near infrared 
and well constrained (Figure]^, and hence we would ex¬ 
pect an extremely high detection significance. We detect 
water at at an extremely high degree of confidence, >17 
cr, in both objects. H 2 S is an example of a poorly con¬ 
strained species (Figure®, and hence would expect, and 
indeed do find, a low detection significance below 2 cr . 
In fact for, HD3651B the Bayes factor is less than one 
(/n(Bayes factor) < 0) or evidence against H 2 S. Since we 
don’t visibly see any obvious spectroscopic features due 
to NH 3 , but do indeed obtain a strong constraint (the 
marginalized posterior is bounded on both sides as op¬ 
posed to an upper limit like H 2 S), we may expect the 
detection significance to fall in between the two afore¬ 
mentioned extremes. This is what we do indeed find, a 
> 6(7 detection of NH 3 in both objects, which is consid¬ 


ered strong. We should note that Bayes factors can be 
sensitive to the prior ranges. We found that in our case, 
the Bayes factor calculation is insensitive to our prior 
ranges. 

We also show two additional, more straight forward, 
lines of evidence in Figure® For the first test we re-ran 
the retrieval but initialized the MCMC with a a non- 
detectable ammonia abundance far from the retrieved 
value. If the retrieved value we obtain is true, then we 
would expect the ensemble of Markov chains to converge 
towards the true value regardless of the starting point- 
given a long enough run time. That is indeed what we 
find. The left two panels in Figure show the evolution 
of the Markov chains. They readily rebound from the 
poor initial starting point (essentially no ammonia) and 
converge to a nice tightly packed bundle within the target 
distribution. Less than 10% of the chains remain outliers. 
Had we run for even longer (weeks perhaps) these chains 
would likely have fallen in line with the rest and too 
converge to the true answer. 

For the final test for ammonia we create two synthetic 
brown dwarf spectra for which we know the true TP pro¬ 
file and composition. We choose parameter values and 
data properties similar to G1570D. We create two syn¬ 
thetic spectra: one with ammonia and one without. We 
then apply the retrieval to these two spectra. What is 
shown in the right panel Figure are the retrieved am¬ 
monia distributions for these two scenarios. The blue 
histogram is the retrieved ammonia distribution for the 
synthetic spectra generated without ammonia. In this 
case only an upper limit of ammonia can be obtained. 
This means that there are no spectral features present in 
the synthetic spectrum to suggest the existence of am¬ 
monia. When the ammonia abundance creeps up to a 
high enough values (in this case ^ Ippm) it begins to 
present itself in the spectrum in an undesirable fashion. 
The red histogram is the retrieved ammonia distribu¬ 
tion for the synthetic spectra generated with ammonia. 
This histogram is nicely bounded on both sides suggest- 
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TABLE 3 

Bayesian nested model comparison supporting the presence 

OF AMMONIA IN THE NEAR IR. THREE MODEL SCENARIOS ARE 
SHOWN, ONE THAT INCLUDES ALL PARAMETERS EXCEPT NH 3 ( 
NH 3 ), ONE THAT INCLUDES ALL PARAMETERS EXCEPT H 2 S (H 2 *S), 
AND ONE THAT INCLUDES ALL PARAMETERS BUT H 2 O (H 2 O). THE 

Bayes factors are computed relative to the full model 

WHICH INCLUDES ALL OF THE PARAMETERS SHOWN IN TABLE ^ . 


Scenario 

G1570D 

InH 

Det. Sig. 

HD365IB 

InH 

Det. Sig. 

NH3 

20.5 

6.7 

20.7 

6.8 

H 2 S 

0.7 

1.8 

-1.7 

- 

H 2 O 

I 66 .I 

18.4 

153.4 

17.7 


ing that the spectral features due to ammonia are enough 
to place both a lower and upper bound on the retrieved 
abundance. These three lines of evidence strongly sug¬ 
gest the presence of ammonia in the the low-resolution 
near infrared spectra of these two T-dwarfs. 

3.2. Verifying Cloud Free 

T-dwarfs are typically assumed to be cloud free given 
their blue colors (Burrows et al. 1997; Allard et al. 2001) 
and the success of cloud-free grid models to reasonably 
explain their spectra (e.g., Stephens et al. 2009), though 
some of the redder objects are better matched with mod¬ 
els that include sulfide-like clouds (Morley et al. 2012). 
Since we assumed cloud free atmospheres for our retrieval 
conclusions, we need to make sure that this is indeed the 
case. We are not interested in the cloud properties them¬ 
selves, rather the impact that some unaccounted for gray 
absorber may have on the spectra. Therefore we model 
clouds rather simplistically as a gray absorber with opac¬ 
ity Kc between two defined pressure levels. Pc,bottom and 
Pc,top- The cloud optical depth is 

Pc,bottom Pc,top /^\ 

Tc = l^c - ( 6 ) 

9 

We assume Pc,bottom — Pc,top spans one scale height so 
that we only need to retrieve the cloud base location 
(Pc,bottom) and Kc. The gray approximation can be read¬ 
ily justified. From grid model investigations with clouds 
(e.g., Morley et al. 2012) high sedimentation values 
(Ackerman & Mar ley 2001) are required to best match 
the spectra. High sedimentation values generally result 
in a wider range of particle sizes thus washing out Mie 
scattering features resulting in gray, or at least, nearly 
gray absorption. 

In order to determine whether or not clouds are 
present, we undergo the sa me B ayesian hypothesis test¬ 
ing procedure described in §3.1[ This time the full model 
includes all of the parameters in Table plus the two 
cloud parameters. We compare the evidence of original 
cloud free model to this new full model. Table [4] shows 
the Bayes factor and detection significance for the cloud. 
We find that for G1570D the detection significance is be¬ 
low 2 a suggesting a weak detection of clouds. HD3651B 
presents a slightly higher, or weak to moderate detection 
of a cloud. Figure shows the parameter distributions 
for both the cloudy and clear atmospheres; including the 
clouds has an insignificant impact on all of the retrieved 
parameter values (e.g., the change in the median value is 
less than the typical width of the distributions). There¬ 
fore, we are justified in assuming cloud free atmospheres 


for these two objects. 


TABLE 4 

Bayesian nested model comparison demonstrating lack of 

EVIDENCE FOR CLOUDS. ThE BAYES FACTORS ARE COMPUTED 
RELATIVE TO THE MODEL WHICH INCLUDES ALL OF THE 
PARAMETERS SHOWN IN TABLE [I] and THE TWO ADDITIONAL CLOUD 
PARAMETERS (SEE TEXT). 



G1570D 


HD365IB 


Scenario 

InB 

Det. Sig. 

InH 

Det. Sig. 

Cloud 

0.76 

1.87 

1.65 

2.38 


3.3. Impaet of Photometry 

The fiux calibrated spectra depend on the choice of 
photometry used as demonstrated in Figure Here, 
we explore the impact on the retrieved quantities of the 
choice in photometry. For each object we calibrate the 
spectrum with either the J-band, H-band, or K-band 
2MASS photometry, as shown in Figure [l] As men¬ 
tioned earlier, Stephens & Leggett (2004) provide cor¬ 
rection factors for the 2MASS photometry. We do not 
apply those correction factors here, rather the goal is to 
determine what effects “bad” photometry may have on 
the retrieved quantities. We then execute the retrieval on 
each of those three calibrated spectra. Figureshows the 
resulting retrieved quantities for each of the photometric 
calibration scenarios. The impact is minimal for most 
quantities (e.g., the shift in the median is well within the 
Icr uncertainties) with the exception of the photometric 
radius. This is unsurprising as the overall scaling to the 
spectrum depends on {R/D)‘^. Shifts in this scaling due 
to photometry will result in changes in the derived ra¬ 
dius. We also find small (^ la) shifts in the retrieved 
gravity due to the prior upper limit on the mass (masses 
cannot exceed 80 Mj). This is because the mass depends 
on both the radius and gravity therefore the small shifts 
in the radius propagate through the mass upper limit 
to the gravity. If the radius increases due to a change 
in photometry, then the gravity has to decrease. These 
shifts in derived radius and gravity stress the importance 
of precision photometry on these objects. 

4. VALIDATING RETRIEVAL WITH BENCHMARKS 

Both G1570D and HD3651B happen to orbit stars with 
known properties. This makes them powerful benchmark 
systems for which we can test the validity of our retrieval 
approach. 

The basic properties we use to evaluate the benchmark 
systems are their evolutionary derived ages, metallici- 
ties, and the carbon-to-oxygen ratios. A fundamental 
assumption with benchmark systems is that the primary 
and companion both formed out of the same nebular ma¬ 
terial at the same time and should each individually in¬ 
dicate the same elemental abundances and age. This is 
unlike planetary systems in which we believe that planet 
formation processes within the protoplanetary disks can 
alter the planetary atmosphere abundances relative to 
their host star (Oberg et al. 2011; Fortney et al. 2013). 

A summary of the relevant stellar primary and re¬ 
trieved brown dwarf properties are shown in Table ^ 
Since stellar photospheres are generally assumed to be 
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Fig. 6 .— Evidence for ammonia in the near infrared. The left two panels show show the evolution of the MCMC chains initialized 
at a non-detectable value. The chains for each object readily converge towards a well constrained solution about the quoted retrieved 
values. The outlier chains account for less than 10% of the total probability. The right panel shows the retrieved probabilitv distribution of 
ammonia for two synthetic brown dwarf spectra. The first synthetic spectrum was generated with a mixing ratio of 10“^®. The retrieved 
probability (blue) shows only an upper limit-this is consistent with a non-detection-as we would expect. The second synthetic spectrum 
was generated with an ammonia mixing ratio similar to the retrieved values for G1570D and HD3651B (vertical dashed line). The retrieved 
probability distribution (red) is bounded on both sides about the truth suggesting a strong constraint. 
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Fig. 7.— Impact of a gray cloud on the retrieved quantities for G1570D (top row) and HD3651B (bottom row). In green we show the 
retrieved marginalized posterior distributions for the cloud free nominal model as in Figure ^and in blue for the cloudy model. The median 
and Ifj confidence interval for each scenario are shown above each histogram. The shift in the medians of all parameters remain less than 
the one sigma uncertainty suggesting that the inclusion of a gray cloud has a minimal impact. 


well mixed and free from condensates, their photo- 
spheric abundances are representative of their intrinsic 
values, however for the brown dwarfs we retrieve the 
atmospheric quantities rather than the intrinsic quan¬ 
tities, which can be different due to the aforemen¬ 
tioned reasons. Therefore, in Table we show both 
the atmospheric elemental quantities and the thermo- 
chemically self-consistent intrinsic elemental quantities 
derived from the hand tuned fits (Figurej^) of the chem¬ 
ical model to the retrieved abundances. Typically the 
atmospheric oxygen is depleted by 20-30% relative to the 
intrinsic due to the sequestration of oxygen in conden¬ 
sates (depending on the intrinsic metallicity and C/0). 
This results in a higher atmospheric carbon to oxygen 
ratio and an overall lower atmospheric metallicity, since 


oxygen is the dominant metal atom. In Table we list 
the inferred intrinsic metallicity and C/0 based on the 
atmospheric measurements and the estimated depletion 
of O and silicates due to condensation. These are the 
nominal values against which we compare to the stellar 
abundances. We find that the metallicities in both ob¬ 
jects are somewhat lower than what is measured in the 
primaries but are still consistent. From the C/0 mea¬ 
surements in both the brown dwarf (Q and host star, 
we have an additional benchmark constraint. We find 
that for the G1570 system, the la inferred intrinsic C/0 
range falls entirely within the stellar primary C/0 val¬ 
ues. We consider this an excellent agreement. For the 
HD3651 system, the inferred intrinsic C/0 is higher than 
the stellar primary by ^30%, however their medians are 
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Fig. 8 . — Impact of choice of photometry used to calibrate the normalized Spex spectra. For each object (G1570D-top row, HD3651B- 
bottom row) we show the marginalized posteriors resulting from the calibrated spectra from each of the three photometric band (blue for 
J-band, green for H-band (nominal case), and red for K-band). We also show the median and Icr confidence interval for each parameter 
for each photometric case. 


consistent at the 2cr level. This represents, for the first 
time, a comparison of a stellar and companion brown 
dwarf carbon-to-oxygen ratios. This suggests that we 
are in fact retrieving the proper molecular abundances 
in the brown dwarf atmospheres. 

Finally, we compare the evolution-derived age to the 
estimated system age. Figure shows the Saumon & 
Marley (2008) isochrones in log ^-effective temperature 
space. The shaded regions represented the estimated 
system ages as described extensively in Liu, Leggett, & 
Chiu (2007). Our retrieved gravity and effective temper¬ 
ature are consistent with the presumed system ages. We 
also obtain photometric masses from the retrieved grav¬ 
ity and photometric radii. We find that the la range in 
photometric masses for G1570D (15-58 M D is consistent 
with the evolution model masses (Figure^; however for 
HD3651B we find somewhat higher photometric masses 
(45-78 Mj) than anticipated from the evolution models. 
It is unclear why this may be. 

In summary, we find ages derived from our retrieved 
gravity and effective temperatures are consistent with the 
measured system age and that the retrieved metallicities 
and C to O ratios, after taking into account the loss of 
oxygen due to condensates, are in good agreement with 
the primaries. 

5. SUMMARY & CONCLUSIONS 

We have established a new minimal-assumption re¬ 
trieval approach for brown dwarf atmospheres that signif¬ 
icantly advances the work of our previous retrieval study 
(Line et al. 2014b). The new retrieval approach relies 
upon a more robust Bayesian estimator, forward model, 
temperature profile parameterization, a single continuous 
spectrum, and treatment of unknown systematic uncer¬ 
tainties permitting generous uncertainty estimates. Dif¬ 
ferences in our results compared with Line et al. (2014b) 
are likely due to these major changes. From our new ap¬ 
proach applied to two benchmark late T-dwarfs, G1570D 
and HD3651B, we determined the allowed range of the 
thermal structures, molecular abundances, gravities, and 
radius-to-distance scalings directly from the data. We 
found that this parameter set provides very good fits to 
the data in the form of minimal, nearly random, resid¬ 
uals. We validated the chemical plausibility of the re- 


TABLE 5 

Benchmark system properties (parameters from Liu, 
Leggett, & Chiu (2007) unless otherwise noted). The C to 
O RATIO FOR Cl570A is from our stellar abundance analysis 
DESCRIBED IN THE APPENDIX. QUANTITIES LABELED WITH 

“retrieved” are the retrieved values from this study. 


Property 

HD365IB 

C1570D 

Spectral Type . 

T7.5 

T7.5 

Host star spectral type . 

KOV 

K4V 

Distance (pc|^ . 

II.06T0.03 

5.84+0.03 

Estimated age (Cyr) . 

3 - 12 

1 - 5 

Host star [Fe/H] . 

O.II - 0.2^ 

-0.22 - 0.1^ 

Retrieved Atmospheric [Fe/H]. 

-0.09 - 0.05 

-0.37 - -0.12 

Chemically Derived Bulk [Fe/Hp] . 

+0.08 

-0.15 

Inferred Intrinsic [Fe/Hp] . 

-0.01 - 0.13 

-0.29 - -0.04 

Host star C/c|^ . 

0.51 - 0.73 

0.65 - 0.97 

Retrieved Atmospheric C/O . 

1.06 - 1.35 

0.95 - 1.25 

Chemically Derived Bulk C/O . 

0.80 

0.70 

Inferred Intrinsic C/c|^ . 

0.76 - 0.97 

0.70 - 0.93 

Retrieved Tg// [K] . 

726^2? 

7141^23 

Retrieved log^f [cm s“^] . 

c; 1 9 +0.09 
-0.17 

4-761°:^^ 

Retrieved Mass [Mj] . 

66 

DO _2i 

31 

-16 

Retrieved Radius [Rj] . 

1 09 

0.08 

1 14 + 0.10 

-0.09 


^van Leeuwen et al. 2007 
^Ramierez et al. (2013) 

^derived from our stellar abundance analysis described in the Ap¬ 
pendix 

•^the chemically derived bulk quantities are the “by hand” ther¬ 
mochemical model fits to the retrieved molecular abundances 
® Assumes an atmospheric metal depletion due to loss of O and 
silicates of 20 % based on the chemical models 
^derived from our stellar abundance analysis described in the Ap¬ 
pendix 

^assuming an O depletion from silicates of 28% for HD365IB and 
26% for C1570D based on the chemical model results 

trieved molecular abundances using a a well vetted ther¬ 
mochemical equilibrium model. 

Perhaps the most significant highlight of our work is 
the robust detection of ammonia in the low resolution 
near-infrared spectra in these late T-dwarfs. We pre¬ 
sented three lines of evidence to support this claim. Fur¬ 
thermore, we showed that clouds play a minimal role in 
sculpting the spectra of these two objects and their inclu¬ 
sion had little to no influence on the other parameters. 
We also suggested that large systematic uncertainties in 
photometry can result in biased estimates of the photo- 
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Fig. 9.— Comparison of our retrieved gravity and effective tem¬ 
perature to the evolution tracks of Saumon & Marley (2008). Our 
retrieved values are the red (G1570D) and blue (HD3651B) boxes 
with error bars. The dotted lines are the logg-T^ff isochrones. 
The red and blue shaded regions are the range of estimated ages 
of the G1570 and HD3651B systems, respectively. The inferred 
evolutionary ages are consistent with the estimated system ages. 


metric radii. 

An additional highlight, is for the first time, using the 
carbon-to-oxygen ratio of the host-companion as an extra 
dimension in establishing benchmark systems. We found 
a remarkable agreement of the carbon-to-oxygen ratios 
derived from a stellar abundance analysis for G1570A 
and our retrieval analysis on G1570D, and a consistent 
agreement within the HD3651 system (considering the 
spread in literature G/0 values). This is quite the ac¬ 
complishment as two completely separate techniques on 
two different objects for which we would expect similar 
abundances, are in good agreement. This further bol¬ 
sters the suggestion of Fortney (2012) to explore the role 
of G/0 in T-dwarf atmospheres and that the G to O ra¬ 
tio should be considered as an additional dimension when 


interpreting brown dwarf spectra. Finally, we found that 
the ages derived from the evolution models and our re¬ 
trieved gravity and effective temperatures are consistent 
with the estimated system ages. It would be interesting 
in future investigations to identify systems for which the 
companion-primary G/0 and metallicities differ by a sig¬ 
nificant amount. This could point to new physics and/or 
chemistry operating in substellar atmospheres. 

This investigation further establishes the power of 
our novel retrieval approach in understanding the atmo¬ 
spheric and bulk properties of brown dwarfs. In a future 
study we plan to apply this technique to a wider range 
of objects with the goal of identifying trends in the ther¬ 
mal structures and molecular abundances and how they 
correlate with empirical metrics. This will undoubt ably 
verify hypothesized physical and chemical mechanisms 
operating in brown dwarf atmospheres, and likely iden¬ 
tify unknown ones as well. 
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7. APPENDIX 

Comparing the carbon-to-oxygen ratio in stellar-brown dwarf companion systems provides an additional benchmark 
dimension. Both brown dwarfs and the stars that they orbit are presumed to form out of the same molecular cloud, 
and thus would be expected to have the same elemental abundances. Metallicity and age are usually the benchmark 
dimensions. Additional dimensions provide more constraints on the system properties. Determining stellar abundances 
is no easy task as different groups using different techniques with different data on the same objects often times report 
significantly different results (Hinkle et al. 2014). In this section we provide our own analysis to determine the stellar 
carbon-to-oxygen ratios in G1570A and HD3651A. 


7.1. Observations and Stellar Parameter Analysis 
7.1.1. G1570A 

The observations of the K4 dwarf G1570A were conducted on 13 July 2014 (UT) with the Magellan Inamori Kyocera 
Echelle (MIKE) spectrograph (Bernstein et al. 2003) on the 6.5m Landon Glay (Magellan II) Telescope at Las 
Gampanas Observatory. Three frames of 100s each were taken of the target with the 0.5”x5” slit and 1x1 binning. 
On 12 July 2014 (UT), three frames of 500s each with the 0.35”x5” slit and 1x1 binning were taken of Vesta, as a 
solar standard. On both nights calibrations (biases, quartz and milky flats, and ThAr lamp spectra) were taken at the 
beginning of the night. MIKE is a double echelle spectrograph, meaning a dichroic splits the light into blue (3350-5000 
A ) and red (5000-9400 A) arms. The data were reduced, extracted, combined, and wavelength calibrated with the 
Garnegie Python Distribution (GarPy) MIKE pipeline, written by D. Kelson (see also Kelson 2003). The resulting S/N 
in the G1570A spectrum was ^160, and ^200 at 6300 A . Gontinuum normalization, order stitching, and Doppler-shift 
correction were performed with standard packages in IRApj^ 

Though the stellar parameters of G1570A have been measured previously (e.g., Eeltzing & Gustafsson 1998; Thoren 
& Eeltzing 2000; Ghezzi et al. 2010; Lee et al. 2011), for consistency we re-derived the Teff, log g, microturbulence 
(^), and [Ee/H] values from the MIKE spectra, based on the methods from our previous work (e.g., Teske et al. 2014). 
Briefly, G1570A’s stellar parameters were derived from equivalent width (EW) measurements of Ee I and Ee 11. We 
used the iron line list of Tsantaki et al. (2013), optimized for cool stars (Teff < 5000 K) by matching spectroscopic and 
infrared flux method (IREM) temperatures. We forced zero correlation between [Ee I/H] and lower excitation potential 


^ IRAF is distributed by the National Optical Astronomy Ob¬ 
servatory, which is operated by the Association of Universities for 


Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 
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(x) to set Teff, zero correlation between [Fe I/H] and reduced equivalent width [log(EW/A)] to set and zero difference 
(within two decimal places) between [Fe I/H] and [Fe II/H] to set log g. The abundances of Fe were determined using 
the local thermodynamic equilibrium (LTE) spectral analysis code MOOG (Sneden 1973), with model atmospheres 
interpolated from the Kurucz ATLAS9 NOVER gridffl Equivalent widths were measured in IRAE with the ‘splot’ 
task, and abundances were normalized to the solar values as measured in our Vesta spectrum on a line-by-line basis. 
The logV(Ee) values for the Sun were determined with our Vesta spectrum and a solar Kurucz model with Teff=5777 
K, log ^=4.44 dex, [Ee/H]=0.00 dex, and ^=1.38 km s“^. In G1570A, 89 Ee I and 10 Ee II lines were measured; the 
line properties and EWs are provided in Table 


TABLE 6 

Lines Measured, Equivalent Widths, and Abundances 


Ion 

A 

(A) 

X 

(eV) 

log gf 

EWq 

(mA) 

logAT© 

EWgL570A 

(mA) 

logA^oLsroA 

EWhD3651A 

(mA) 

logA^HDsesiA 

C I 

5052.17 

7.685 

-1.24 

33.3 

8.38 



25.0 

8.75 

C I 

5380.34 

7.695 

-1.57 

20.9 

8.44 



15.1 

8.81 

[C I] 

8727.13 

1.26 

-8.165 

5.3 

8.43 

6.8 

8.568 



[O I] 

6300.30 

0.00 

-9.717 

5.4 (5.0) 

8.67H 

8.6p 

(8.62“) 

b-fil*’) 

7.9 

8.59“ 

8.54^ 

6.6 

8.82“ 

8.82^ 

O I 

7771.94 

9.15 

0.369 

71.5 (69.7) 

8.8^ 

(8.83=) 

12.3 

8.85“ 

38.4 

9.14“ 


7774.17 

9.15 

0.220 

61.5 (63.2) 

8.86^ 

(8.88=) 

11.2 

8.93“ 

38.3 

9.28“ 


7775.39 

9.15 

0.001 

49.0 54.8) 

8.86^ 

(8.78=) 

7.30 

8.87“ 

25.1 

9.15“ 


Note. — The number abundances {\ogN) for HD3651A listed in this table are calculated as an example with Allende Prieto et al. 
(2004) stellar parameters. Any value in parentheses refers to a HIRES solar measurement; all solar-normalized HD3651A abundances in 
the text are relative to HIRES solar measurements. This table is available in its entirety in a machine-readable form in the online journal. 
A portion is shown here for guidance regarding its form and content. 

^Abundance derived through equivalent width analysis. 

^Abundance derived through synthesis analysis. 

LTE abundance. 


Our final parameters and errors for G1570A are listed below (Table along with those of several other studies for 
comparison. The errors are calculated as in our previous work (Teske et al. 2013ab, 2014, 2015) - the change in 
Teff (^) required to cause a correlation coefficient r between [Ee I/H] and y ([Ee I/H] and reduced EW) significant at 
the Icr level was adopted as the uncertainty in these parameters. The uncertainty in log g was calculated differently, 
through an iterative process described in detail in Baubar & King (2010). Uncertainties in [Ee I/H] and [Ee H/H] are 
calculated from the quadratic sum of the individual uncertainties in these abundances due to the derived uncertainties 
in Teff, log g^ and as well as the uncertainty in the mean (cr^ of each abundance. The uncertainty due to the stellar 
parameters is measured from the sensitivity of the abundanc^o each parameter for changes of ±150 K in Teff, ±0.25 
dex in log g, and ±0.30 km s“^ in The uncertainty due to each parameter is then the product of this sensitivity and 
the corresponding parameter uncertainty. Eor the abundances determined through spectral synthesis (e.g., from [O 
I], see below), models with this range of stellar parameters were compared to the data and the elemental abundance 
adjusted to determine the best fit. 

While our log g value is moderately lower than other some other studies, it agrees within errors. As described below, 
these stellar parameter errors are propogated through the other abundance measurements. 


TABLE 7 

Gl570A Stellar Parameters 


Parameter 

this work 

Eeltzing 

& Gustafsson (98) 

Thoren 

& Eeltzing (00) 

Ghezzi et al. (10) 

Lee et al. (11) 

Teff (K) 

4686 ±47 

4585 

4585 

4799 ±72 

4615 

log g (cgs) 

4.37 ±0.27 

4.70 

4.58 

4.60 ±0.16 

4.36 

^ (km s“^) 

1.03 ±0.16 

1.0 

1.0 

0.77±0.08 


[Fe/H] (dex) 

-0.05±0.17 

0.04 

0.04 

0.03±0.03 

-0.05 


7.1.2. HD3651A 

The bright (U=5.88) KO dwarf HD3651A hosts a 0.2 Mj planet (Eischer et al. 2003), and has thus been the 
target of many spectroscopic observations and stellar parameter analyses (^15, according to SIMBAD). Given the 
many previous stellar parameter analyses based on high-resolution, high-S/N data, we do not derive yet another set 
of stellar parameters here. Instead we use an archive HIRES spectrum (J. Johnson 2015, private communication). 


^ See http://kurucz.harvard.edu/grids.html 

^ o-fj, = cr/VA/”—1, where cr is the standard deviation of the 
derived abundances and N is the number of lines used to derive 


the abundance. 
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with S/N^200 at the [O I] 6300 A line, along with several reported sets of stellar parameters (Table to verify the 
previously-measured carbon and oxygen abundances. We also assess realistic uncertainties on these measurements, as 
no formal uncertainties were previously published; our analysis is reported in the next section. The solar standard in this 
case is an archive HIRES spectrum of reflected light from the asteroid Vesta (A. Howard 2014, private communication), 
taken in the same configuration as the HD3651A spectrum. 

TABLE 8 

HD3651A Previously Measured Stellar Parameters 


Parameter 

Allende Prieto et al. (04) 

Delgado Mena et al. (10) 

Petigura & Marcy (11) 

Ramirez et al. (13) 

Teff (K) 

5117A94 

5173 

5221 

5303±63 

log g (cgs) 

4.58 ±0.04 

4.37 

4.45 

4.56±0.03 

i (km s-l) 

0.93 

0.74 


0.64±0.12 

[Fe/H] (dex) 

0.13 

0.12 

0.16 

0.18±0.07 


7 . 2 . C/0 Ratio Analysis 
7 . 2 . 1 . Carbon Measurements 
7 . 2 . 2 . G1570A 

The determination of [C/H] in G1570A required a different technique than our previous work, due to the low 
temperature of the star. The carbon abundances for G1570A derived from EW measurements of widely-used high- 
excitation (x > 7.67 eV) G I lines (5380, 7711, 7113A ) and a line-by-line comparison to solar resulted in [G/H]avg=0.99, 
an unrealistically high value. These G I lines suffer NLTE effects, such that an LTE analysis overestimates the 
abundances, and corrections based on NLTE atomic models are predicted to increase in magnitude (larger negative 
corrections) with higher Teff and lower log g (e.g., Asplund 2005; Takedafe Honda 2005; Eabbian et al. 2006). However, 
similar to the O I triplet lines at 7771/7774/7775A (see discussion below), the NLTE corrections for cool stars are 
predicted to be minimal, less than the solar-type star corrections of ^-0.05. As discussed in Teske et al. (2013a) 
for the cool star (^5350 K) 55 Gnc, applying the predicted NLTE corrections to the O I triplet abundances actually 
increases the [0/H] values, rather than decreases, because the solar corrections exceed the cooler-star corrections and 
thus their differences are increased. We find the same problem with the G I lines in G1570A. 

Thus, the [G/H] for G1570A is instead derived from the low-excitation (x=1.26 eV) forbidden [G I] line at 8727.13A 
(Eigure[lQ|. This line is weak, and possibly blended with an Ee I line at 8727.lOA (Lambert & Swing 1967), but has a 
well-determined transition probability and is not susceptible to depatures from LTE (Gustafsson et al. 1999; Asplund 
et al. 2005). Using a spectral synthesis analysis with a line list between 8724-8730A , gathered from the Vienna Atomic 
Line Database (VALD; Kupka et al. 1999), we derive A(G)gl 570 A= 8 . 57 . This also agrees with our derivation using 
the IRAE-measured EW (5.30 mA ) and the ‘abfind’ driver in MOOG with the G1570A model derived from the same 
MIKE spectra. Via the same procedure, we measure A(G)soiar=8.43 from the 8727 A [G I] line. The resulting [G/H] is 
listed in TableNote that in Tablewe include the formal errors on each abundance measurement ([Ni/H], [G/H], 
[0/H]), which IS the quadratic sum oithe three individual parameter uncertainties (Teff, log f) and cr^ as described 
above for the case of iron. 


7 . 2 . 3 . HD3651A 

The [G I] line used to measure the carbon abundance of G1570A was not available in our HD3651A HIRES spectrum, 
so we instead rely on the lowest excitation (x=7.685 eV) available G I lines at 5052.2 A and 5380.3 A. These lines, as 
mentioned above, are predicted to have negligible NLTE corrections at the low temperature of HD3651A (^5100-5200 
K), on par with those predicted for the Sun (<0.05 dex; Takeda & Honda 2005). Thus, any deviations from LTE 
should cancel with the calculation of the solar-normalized [G/H] value derived from these lines. In this case, the G 
I lines resulted in reasonable [G/H] values (see Tables and [1Q|) . Measurements of other available G I lines used in 
previous work (6588, 7111, 7113A, e.g. Teske et al. 2014) result in ^0.1-0.25 dex larger [G/H] values - as expected 
since these higher excitation lines are more susceptible to NLTE effects - and thus are excluded for the final [G/H] for 
HD3651A. We use our equivalent width measurements (Table |6) with three different stellar models (Allende Prieto et 
al.’s, Delgado Mena et al.’s, and Petigura & Marcy’s, see Table ^ to derive carbon abundances, and use the resulting 
spread in [G/H] as a measure of its uncertainty. Our measureo/G/H] values (0.21-0.37) agree with those of Allende 


Prieto et al. (2004) and Delgado Mena et al. (2010), as expected (Table p^. Note that in Table [TOj we include only 
the standard deviation (a) errors, since the stellar parameter errors are not derived in this work, and are calculated 
differently for each source; the error calculation method used for G1570A does not apply. 


7 . 2 . 4 . Oxygen Measurements 

As discussed in Teske et al. (2013) and (2014), and by many others in the past (e.g., Nissen & Edvardsson 1992; 
King & Boesgaard 1995; Asplund et al. 2004; Schuler et al. 2006a,b; Gaffau et al. 2008), oxygen abundances are 
notoriously difficult to derive, particularly in stars that are cooler/hotter and/or more/less metal-rich than the Sun. 
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The oxygen abundance indicators we explored here include the forbidden [O I] line at 6300A, which is well-described 
by LTE but blended with a Ni I line, and the O I triplet at 7771-7775A , made up of three unblended and usually 
prominent lines amenable to direct EW measurement. 

We derived the [O I] 6300.304A line using two methods. Eirst, we measured the EW of the feature directly from 
the spectrum as input for the ‘blends’ driver in MOOG, accounting for ^^Ni I + ^^Ni I feature at 6300.335A with 
log5f/(^^Ni)=-2.695 and log5f/(^^Ni)=-2.275 as derived by Bensby et al. (2004). This resulted in [O/Hjesoo,blends = -0.08 
for G1570A and 0.19-0.23 for HD3651A, depending on the set of stellar parameters. Second, we performed a spectral 
synthesis analysis on the line, using as a “known” our measured [Ni/H] abundance based on a line-by-line analysis with 
the Sun (as with Ee; see Tablesj^and 10). The synthesized spectra (with the stellar parameters derived or noted above) 
were convolved with a Gaussian prome, based on near-by unblended lines, to represent the instrument PSE, stellar 
macroturbulence, and rotational broadening; we also fixed the nickel abundance to our measured value. Unlike in the 
case of 55 Gnc (Teske et al. 2013a), the oxygen line strength did not change drastically (^0.02 dex) by changing the 
Ni abundance within our derived [Ni/H]. The remaining free parameters were continuum normalization, wavelength 
shift, and oxygen abundance. The best fit to the synthesized spectra for the [O I] was determined by minimizing the 
deviations between the observed and synthetic spectra (see Eigure 10). This resulted in [ 0 /H] 63 oo,synth = -0.13 for 
G1570A and 0.17-0.27 for HD3651A, depending on the set of stellar parameters. Eor our final [O/Hjesoo we took the 
mean of these measurements, -0.11T0.19 dex for G1570A and 0.18-0.25 dex for HD3651A. 

The O I triplet suffers NLTE effects due to the dilution of the each line’s source function with respect to the 
Planck function (e.g., Kiselman 1993; Gratton et al. 1999; Kiselman 2001), so abundances derived assuming LTE are 
overestimated. As with G I, the predicted NLTE corrections increase with decreasing gas pressures and/or increasing 
temperatures, but decrease for cool stars (e.g., Takeda 2003; Ramirez et al. 2007, Eabbian et al. 2009) like both G1570A 
and HD3651A. We measured the EWs of the three O I triplet lines and derived an A(0) for each line in G1570A, 
HD3651A, and their respective solar standards (see Table 1^; the average [0/H]lte of the line-by-line differences with 
the Sun is 0.05±13 for G1570A (line 3 of Table ||) and 7-0.38 for HD3651 (lines 3, 6, and 9 of Tabled. We 
then applied NLTE corrections from three different sources - Takeda (2003), Ramfrez et al. (2007), and EaBbian et 
al. (2009) - to G1570A, HD3651, and the respective solar standard measurements, and recalculated the line-by-line 
abundance differences (see discussion in Teske et al. 2013a for details of each of these correction schemes). Gorrections 
in absolute abundance (not relative to solar) for the Sun range from 0.13-0.21, for G1570A range from 0.04-0.07, and 
for HD3651 range from 0.52-0.85. The resulting [O/ILnlte, averaged over all three lines and all three sources of NLTE 
corrections, is 0.14 dex for G1570A (line 4 of Table and 0.25-0.46 dex, depending on the stellar parameters. In the 
case of G1570A, based on previous NLTE abundance uncertainty calculations, we assume the same uncertainty as the 
LTE abundance; for HD3651, again the a (standard deviations) across the three lines and three sources are listed in 
Table [Tol 


TABLE 9 

Gl570A Abundances & Indicators with Formal Errors 


Source 

[Ni/H] (dex) 

[C/H] (dex) 

[O/H] (dex) 

[O/H] indicator 

C/O 

this work 

0.01 ±0.05 
0.01 ±0.05 
0.01 ±0.05 

0.14±0.11 

0.14±0.11 

0.14±0.11 

-0.11±0.19 

0.05±0.13 

0.14±0.13 

[O I] 6300 A 

O I triplet 7775 A LTE 

O I triplet 7775 A NLTE 

0.97±0.22 

0.68±0.17 

0.55±0.17 

Feltzing & Gustafsson (98) 

0.16 

0.18 ±0.18 


[O I] 6300 A 


Petigura & Marcy (11) 

0.15 ±0.06 

0.18 


[O I] 6300 A 



TABLE 10 

HD3651A Abundances & Indicators with ct Errors 


Source 

[Ni/H] (dex) 

[C/H] (dex) 

[O/H] (dex) 

[O/H] indicator 

G/O 

this work, 

0.19±0.04 

0.37±0.002 

0.21±0.01 

[O I] 6300 A 

0.79 

AP04 Params 

0.19±0.04 

0.37±0.002 

0.36±0.05 

O I triplet 7775 A LTE 

0.56 


0.19±0.04 

0.37±0.002 

0.44±0.05 

O I triplet 7775 A NLTE 

0.47 

this work, 

0.22 ±0.05 

0.25±0.004 

0.18±0.02 

[O I] 6300 A 

0.65 

DM10 Params 

0.22 ±0.05 

0.25±0.004 

0.24±0.05 

O I triplet 7775 A LTE 

0.56 


0.22 ±0.05 

0.25±0.004 

0.31±0.05 

O I triplet 7775 A NLTE 

0.48 

this work. 

0.19 ±0.04 

0.21±0.001 

0.25±0.03 

[O I] 6300 A 

0.50 

PM 11 Params 

0.19 ±0.04 

0.21±0.001 

0.16±0.05 

O I triplet 7775 A LTE 

0.62 


0.19 ±0.04 

0.21±0.001 

0.23±0.04 

O I triplet 7775 A NLTE 

0.52 

Allende Prieto et al. (04) 

0.27 

0.26 

0.23 

[O I] 6300 A 

0.59 

Delgado Mena et al. (10) 

0.15 

0.25 

0.06 

[O I] 6300 A 

0.85 

Petigura & Marcy (11) 

0.24 


0.07±0.08 

[O I] 6300 A 


Ramirez et al. (13) 



0.05± 0.04 

O I triplet 7771-5 A LTE 





0.12±0.04 

O I triplet 7771-5 A NLTE 
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[C I] 8727 Line 



[0 I] 6500 Line 



0 I Triplet 



Fig. 10.— Plotted are the the lines measured in this work to determine the C/O ratio of G1570A. Dots represent the spectrum of G1570A, 
while triangles represent the solar standard spectrum. In the top two plots, synthesis fits to the lines are shown in red dot-dashed lines. 
The result abundances are given in Table 

7.3. What are the C/O Ratios of G1570A and HD3651A? 

We calculate the C/O ratid^f G1570A and HD3651A with the Asplund et ah (2009) solar A(C)=8.43 and A(0)=8.69 
values and our [C/H] and [O/K] values as follows: 

C/O = 

with the error on C/O represented by the errors of [C/H] and [0/H] added in quadrature. From our different [0/H] 
indicators, the C/O ratio for G1570A ranges from 0.55±0.17 to 0.97 ±0.22, and for HD3651A ranges from 0.47-0.79, 
depending on the set of stellar parameters and oxygen abundance indicators. 

The [O I] line has been designated as a consistently reliable oxygen abundance indicator (e.g., Lambert 1978; Allende 
Prieto et al. 2001; Asplund et al. 2004; Schuler et al. 2006a) due to its formation in the ground state making an 
LTE approximation exceedingly good (Gaffau et al. 2008). As noted above, the O I triplet suffers significant NLTE 
effects, which are predicted to decrease with temperature. However, studies by Schuler et al. (2004,2006b) and King & 
Schuler (2005) of dwarf stars in the Pleiades, M34, and Hyades open clusters, and the Ursa Major moving group found 
that [0/H] triplet,LTE valucs derived from the O I triplet significantly inereased with decreasing Teff (<5400). If the 
assumption holds that stars within a single cluster or moving group should be chemically homogenous, the increasing 
[0/H] triplet,LTE with Teff is in direct contrast with all the available NLTE calculations. These studies do not point 
to a definitive cause of the NLTE correction discrepancy in cool dwarfs, though they suggest age likely plays a role. 
However, the studies do indicate for cool stars like G1570A and HD3651A, the O I triplet-tempearture trend appears 
to contradict the predicted NLTE oxygen abundance corrections. Due to larger corrections for the Sun versus G1570A 
and HD3651A, applying the NLTE-corrections of Takeda (2003), Ramfrez et al. (2007), and Eabbian et al. (2009) all 


® The G/O ratio - the ratio of the number of carbon atoms to 
oxygen atoms - is calculated in stellar abundance analysis as G/0= 


where log(Nx)=logio(Nx/NH)±12. 
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result in [0/H] values for G1570A and HD3651 that are in general larger than in the LTE case, and C/0 values that 
are in general smaller. 

Given the caveats of the O I triplet NLTE values, we chose here to combine the [0/H] triplet ,lte values from the three 
O I lines, and [O/Hjesoo foi* our final best estimate for the oxygen abundances of G1570A and HD3651A, resulting 
in [0/H]avg = -0.03±0.12 for G1570A (where the uncertainty here is the errors of the O I and [O I] abundances 
added in quadrature) and [0/H]avg =0.21-0.29 for HD3651A, depending on the set of stellar parameters. With 
C/Hgl 570 A =0.14±0.11, our final C/0 for G1570A =0.81±0.16. For HD3651, taking [0/H]avg =0.23, a =0.07 and 
C/H]avg =0.28, a =0.08 from the three different stellar parameter analyses, the final C/Oavg for HD3651A=0.62, 
cr =0.11. 


7 . 4 . Comparison with the Literature 
7 . 4 . 1 . G1570A 

Our [0/H]avg differs from the [0/H] values for G1570A reported by Feltzing & Gustafsson (1998) and Petigura & 
Marcy (2011), who both find [0/H] ^0.15. This higher oxygen abundance, combined with our measured [C/H] (none 
of the other studies measured carbon in G1570A), would lower the C/0 ratio to 0.54, which matches the solar value 
(C/Oq= 0.55±0.10; Asplund et al. 2009; Caffau et al. 2011). Considering our errors, and those reported by Petigura & 
Marcy (2011) for [0/H], the high and low C/0 values would just barely overlap within errors. 

Both Feltzing & Gustafsson (1998) and Petigura & Marcy (2011) also measure higher [Ni/H] values (0.18 dex), 
although the Feltzing & Gustafsson abundance overlaps with ours within errors. Using a stellar model with Pe¬ 
tigura & Marcy (2011)’s derived parameters for G1570A (4744 K Teff, 4.76 dex log g, 0.10 dex [Fe/H], and assuming 
^ =1.00 km s“^) and our measured EWs for Ni I, C I, [O I], and the O I triplet in LTE results in [Ni/H]=0.14 
and C/O=10^‘^^+‘^^/10^‘^^+^^‘^^+^‘^^^/^ =0.89. Performing the same exercise with Feltzing & Gustafsson (1998)’s 
derived stellar parameters (4585 K Teff, 4.70 dex log g, 0.04 dex [Fe/H], 1.00 km s“^^) results in [Ni/H]=0.14 and 
C/O=10^‘^^+‘^^/10^‘^^+^^‘^^+^‘^^^/^ =0.63, where in both formulae the average [0/H]=([0/H]tripiet,LTE+[0/H]63oo)/2. 
Given our uncertainties and the typical uncertainties in [C/H] and [0/H], particularly in cool stars, these alternative 
C/0 values and our reported C/0 agree within errors. This exercise also demonstrates that the different [0/H ]5300 
and [Ni/H] abundances derived in this work versus in Feltzing & Gustafsson (1998) and Petigura & Marcy (2011) are 
likely due to differences in stellar parameters, and not the quality of the spectra or the empirical measurements - using 
our measurements and their models results in [0/H ]5300 and [Ni/H] values very similar to what they report. 

7 . 4 . 2 . HD3651A 

The [C/H] values we derive for HD3651A based on our measurements of a high-S/N archive HIRES spectrum are 
in decent agreement with those reported by Allende Prieto et al. (2004) and Delgado Mena et al. (2010). Since our 
[C/H] is derived from stronger C I lines, versus Allende Prieto’s [C I] 8727 A measured [C/H], we might expect our 
abundance to be higher, which it is. We have perfect agreement with Delgado Mena et al. (2010) when using their 
stellar parameters and our C I equivalent width measurements ([C/H] =0.25). 

The most challenging aspect of the C/0 measurement in this (and many) stars is pinning down the [0/H]. In 
comparison to Allende Prieto et al. (2004)’s [0/H] derived from the [O I] 6300 A line, our [O I] 6300 A measurement 
combined with their stellar parameters produces almost an almost identical [0/H] (0.21 versus 0.23). This is not the 
case when comparing [0/H]5300 values of Delgado Mena et al. (2010) and Petigura & Marcy (2011) to those measured 
here, where we find oxygen abundance higher by 0.12 and 0.18, respectively. Other authors (Fortney 2012; Nissen 
2013; Teske et al. 2013a; Nissen et al. 2014) have called into question the high C/0 ratios (often caused by low oxygen 
abundances) reported in the previous papers, and in some cases have reported different [0/H] results for the same 
stars. Differences in [0/H] measured from the same line in the same star could arise due to several challenging aspects 
of the 6300 A line, such as continuum placement, telluric contamination, weakness of the line at low metallicity, and 
blending with an Ni I line that can make up to 30% of the line strength in the Sun (Caffau et al. 2008). Delgado 
Mena et al. (2010) specifically removed spectra with obvious telluric contamination, and estimated the EW of the Ni 
line blended with [O I] to form the 6300 A line using the “ewfind” driver of MOOG and the Ni abundances measured 
from ^50 Ni I lines. The [Ni/H] value we derived (0.22) from 27 Ni I lines and the same stellar parameters at Delgado 
Mena et al. is slightly lower than their value (0.15). Presumably a smaller Ni abundance comes from a smaller Ni EW, 
meaning a smaller contribution to the 6300 A line and thus a larger contribution from O, and yet Delgado Mena et al. 
also find a smaller [0/H]. The Ni I and O I line parameters that we use (Bensby et al. 2004 for Ni I; Story & Zeippen 
2000 for [O I]) differ from those used by Delgado Mena et al. (Allende Prieto et al. 2001, for Ni I; Lambert et al. 1978 
for [O I]), but this is not enough to account for the 0.12 dex difference. Without a published EW measurements or 
synthesis fits, it is unclear why our [0/H]5300 for HD3651A differs from Delgado Mena et al.’s. 

In their [O I] 6300 A measurements, Petigura & Marcy (2011) also discard any stars with telluric line contamination, 
but many of their spectra are contaminated by iodine lines, which are ^5% deep in the region of the line. To account 
for this, they shift the stellar spectrum (with iodine) to match in wavelength the most recent iodine reference spectrum, 
and then divide the stellar spectrum by the iodine reference spectrum. The resulting spectrum can contain artifacts at 
the ^1% level, likely within any EW measurement. Petigura & Marcy treat the Ni blend differently: Once they derive 
an [0/H] from the 6300 A line by comparing synthetic spectra to their observed spectra, they refit the oxygen line to 
spectra with ±0.03 dex Ni, and add the resulting errors in [0/H] in quadrature to their statistical errors. Petigura 
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& Marcy authors also use different line parameters for the Ni I line blended with oxygen at 6300 A), determined by 
fitting their solar spectrum to match their adopted solar abundance distribution. The [Ni/H] abundance derived by 
Petigura & Marcy for HD3651A (0.24) agrees with what we derive using our EWs measured from 27 Ni I lines and their 
stellar parameters, but they do not explicitly use the measured nickel abundance for each star in their measurement 
of [O/Hjesoo- This may be the reason behind our differing [O/Hjesoo values. 

Ramfrez et al. (2013) do not measure [0/H] from the 6300 A forbidden line, but instead from the O I 7771-5 A triplet. 
As noted above, these lines suffer NLTE effects that are not well understood or calibrated for cool stars. Combining 
our triplet line EWs with the stellar parameters of HD3651A from Ramfrez et al. (who do not list their measured 
EWs) results in [O/H]tripiet,LTE=0.11, cr=0.05, which is slightly higher than, but overlaps within errors of, Ramfrez 
et al.’s [0/H]tripiet,LTE=0.05±0.04 (where here we have added Ramfrez et al.’s line-to-line scatter, 0.01, and their 
uncertainty in the stellar parameters, 0.04, in quadrature for a total error). Similarly, we find [O/H]tripiet,NLTE=0.19, 
<7=0.04, using our triplet EWs, Ramfrez et al.’s stellar parameters and Ramfrez et al. (2007) ’s NLTE correction 
scheme, whereas Ramfrez et al. (2013) reports [0/H]tripiet,LTE=0.12±0.04 for HD3651A. Thus, while slightly higher, 
our [0/H] values still consistent with Ramfrez et al.’s within errors; if we use [0/H] triplet ,lte based on Ramfrez et 
al.’s stellar parameters, and recalculate [C/H] with the same parameters, the resulting HD3651 C/O=0.66, in good 
agreement with our final value above (0.62±0.11). 

Facilities: Magellan: Clay 
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